Tarea 1: Foundations

CC6104: Statistical Thinking

Integrantes :

Nicolás Acevedo Rivas
Diego Navarrete Gonzalez

Cuerpo Docente:

Profesor: Felipe Bravo M.
Auxiliares: Martín Paredes, María José Zambrano
Ayudantes: Nicolas Cabello, Kevin Iturra, Mariana Vásquez

Índice:

Objetivo
Instrucciones
Referencias
Segunda Parte: Elaboración de Código

Objetivo

Bienvenid@s a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en análisis exploratorio de datos y conceptos introductorios de probabilidades. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de ejercicios prácticos con el fín de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

La tarea se realiza en grupos de máximo 2 personas. Pero no existe problema si usted desea hacerla de forma individual.
La entrega es a través de u-cursos a más tardar el día estipulado en la misma plataforma. A las tareas atrasadas se les descontará un punto por día.
El formato de entrega es este mismo Rmarkdown y un html con la tarea desarrollada. Por favor compruebe que todas las celdas han sido ejecutadas en el archivo html.
Al momento de la revisión tu código será ejecutado. Por favor verifica que tu entrega no tenga errores de compilación.
No serán revisadas tareas desarrolladas en Python.
Está PROHIBIDO la copia o compartir las respuestas entre integrantes de diferentes grupos.
Pueden realizar consultas de la tarea a través de U-cursos y/o del canal de Discord del curso.

Referencias:

Slides de las clases:

Videos de las clases:

Introduction to Statistical Thinking: video1 video2
Introduction to R: video1 video2 video3 video4
Descriptive Statistics: video1 video2 video3 video4
Probability: video1 video2 video3 video4 video5 video6

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Pregunta 1: Visualización de Datos

Para esta pregunta usted deberá trabajar en base al conjunto de datos hearth_database.csv, el cual esta compuesto por las siguientes variables:

target: Señala si el paciente tuvo un infarto.
sex: Sexo de los sujetos de prueba.
fbs: Azúcar en la sangre con ayunas. Esta variable señala solo si se encuentra <=120 o >120.
exang: Angina de pecho inducida por el ejercicio.
cp: Tipo de dolor de pecho.
restecg: Resultados electrocardiográficos en reposo.
slope: Pendiente del segmento ST máximo de ejercicio.
ca: Número de buques principales.
thal: Thalassemia.
age: Edad en años.
trestbps: Presión arterial en reposo.
chol: colesterol sérico en mg/dl.
thalach: Frecuencia cardíaca máxima alcanzada.
oldpeak: Depresión del ST inducida por el ejercicio en relación con el reposo.

En base al dataset propuesto realice un análisis exploratorio de los datos (EDA). Para su análisis enfoquen el desarrollo en las siguientes tareas:

Obtenga la media, mediana, quintiles y valores máximos desde los datos que componen el dataset.
Obtenga la Matriz de correlación de Pearson y visualice la relación entre las variables numéricas.
Visualice los boxplot para las variables numéricas.
Visualice a través de un histograma como distribuyen las variables respecto a los TARGET.

Respuesta

my.frame <- read.table(file = "hearth_database.csv",header = T,sep = ",", fileEncoding = "latin1")
colnames(my.frame)[colnames(my.frame) == "ï..age"] <- "age"


#slope
mean(my.frame$slope)

## [1] 1.39934

median(my.frame$slope)

## [1] 1

quantile(my.frame$slope)

##   0%  25%  50%  75% 100% 
##    0    1    1    2    2

max(my.frame$slope)

## [1] 2

#ca
mean(my.frame$ca)

## [1] 0.7293729

median(my.frame$ca)

## [1] 0

quantile(my.frame$ca)

##   0%  25%  50%  75% 100% 
##    0    0    0    1    4

max(my.frame$ca)

## [1] 4

#thal
mean(my.frame$thal)

## [1] 2.313531

median(my.frame$thal)

## [1] 2

quantile(my.frame$thal)

##   0%  25%  50%  75% 100% 
##    0    2    2    3    3

max(my.frame$thal)

## [1] 3

#age
mean(my.frame$age)

## [1] 54.36634

median(my.frame$age)

## [1] 55

quantile(my.frame$age)

##   0%  25%  50%  75% 100% 
## 29.0 47.5 55.0 61.0 77.0

max(my.frame$age)

## [1] 77

#trestbps
mean(my.frame$trestbps)

## [1] 131.6238

median(my.frame$trestbps)

## [1] 130

quantile(my.frame$trestbps)

##   0%  25%  50%  75% 100% 
##   94  120  130  140  200

max(my.frame$trestbps)

## [1] 200

#chol
mean(my.frame$chol)

## [1] 246.264

median(my.frame$chol)

## [1] 240

quantile(my.frame$chol)

##    0%   25%   50%   75%  100% 
## 126.0 211.0 240.0 274.5 564.0

max(my.frame$chol)

## [1] 564

#thalach
mean(my.frame$thalach)

## [1] 149.6469

median(my.frame$thalach)

## [1] 153

quantile(my.frame$thalach)

##    0%   25%   50%   75%  100% 
##  71.0 133.5 153.0 166.0 202.0

max(my.frame$thalach)

## [1] 202

#oldpeak
mean(my.frame$oldpeak)

## [1] 1.039604

median(my.frame$oldpeak)

## [1] 0.8

quantile(my.frame$oldpeak)

##   0%  25%  50%  75% 100% 
##  0.0  0.0  0.8  1.6  6.2

max(my.frame$oldpeak)

## [1] 6.2

#summary(my.frame)

cor(my.frame[,7:14])

##                slope          ca        thal         age    trestbps
## slope     1.00000000 -0.08015521 -0.10476379 -0.16881424 -0.12147458
## ca       -0.08015521  1.00000000  0.15183213  0.27632624  0.10138899
## thal     -0.10476379  0.15183213  1.00000000  0.06800138  0.06220989
## age      -0.16881424  0.27632624  0.06800138  1.00000000  0.27935091
## trestbps -0.12147458  0.10138899  0.06220989  0.27935091  1.00000000
## chol     -0.00403777  0.07051093  0.09880299  0.21367796  0.12317421
## thalach   0.38678441 -0.21317693 -0.09643913 -0.39852194 -0.04669773
## oldpeak  -0.57753682  0.22268232  0.21024413  0.21001257  0.19321647
##                  chol      thalach     oldpeak
## slope    -0.004037770  0.386784410 -0.57753682
## ca        0.070510925 -0.213176928  0.22268232
## thal      0.098802993 -0.096439132  0.21024413
## age       0.213677957 -0.398521938  0.21001257
## trestbps  0.123174207 -0.046697728  0.19321647
## chol      1.000000000 -0.009939839  0.05395192
## thalach  -0.009939839  1.000000000 -0.34418695
## oldpeak   0.053951920 -0.344186948  1.00000000

#boxplot(my.frame[,7:14],main="Hearth")
boxplot(my.frame$slope,main="Hearth slope")

boxplot(my.frame$ca,main="Hearth ca")

boxplot(my.frame$thal,main="Hearth thal")

boxplot(my.frame$age,main="Hearth age")

boxplot(my.frame$trestbps,main="Hearth trestbps")

boxplot(my.frame$chol,main="Hearth chol")

boxplot(my.frame$thalach,main="Hearth thalach")

boxplot(my.frame$oldpeak,main="Hearth oldpeak")

hist(my.frame$age[my.frame$target == "YES"],main="Distribución de Edad para 'YES'")

hist(my.frame$age[my.frame$target == "NO"],main="Distribución de Edad para 'NO'")

hist(my.frame$slope[my.frame$target == "YES"],main="Distribución de slope para 'YES'")

hist(my.frame$slope[my.frame$target == "NO"],main="Distribución de slope para 'NO'")

hist(my.frame$ca[my.frame$target == "YES"],main="Distribución de ca para 'YES'")

hist(my.frame$ca[my.frame$target == "NO"],main="Distribución de ca para 'NO'")

hist(my.frame$thal[my.frame$target == "YES"],main="Distribución de thal para 'YES'")

hist(my.frame$thal[my.frame$target == "NO"],main="Distribución de thal para 'NO'")

hist(my.frame$trestbps[my.frame$target == "YES"],main="Distribución de trestbps para 'YES'")

hist(my.frame$trestbps[my.frame$target == "NO"],main="Distribución de trestbps para 'NO'")

hist(my.frame$chol[my.frame$target == "YES"],main="Distribución de chol para 'YES'")

hist(my.frame$chol[my.frame$target == "NO"],main="Distribución de chol para 'NO'")

hist(my.frame$thalach[my.frame$target == "YES"],main="Distribución de thalach para 'YES'")

hist(my.frame$thalach[my.frame$target == "NO"],main="Distribución de thalach para 'NO'")

hist(my.frame$oldpeak[my.frame$target == "YES"],main="Distribución de oldpeak para 'YES'")

hist(my.frame$oldpeak[my.frame$target == "NO"],main="Distribución de oldpeak para 'NO'")

Pregunta 2: Teorema Central del Limite

Pruebe el teorema central del limite aplicando un muestreo de la media en las distribuciones Gamma, Normal y una a su elección. Grafique los resultados obtenidos y señale aproximadamente el numero de muestreos necesarios para obtener el resultado esperado, pruebe esto con las siguientes cantidades de muestreo \(\{10,100,1000,5000\}\). ¿El efecto ocurre con el mismo número de muestreo para todas las distribuciones?.

Respuesta

# Definición de variables o estructuras necesarias para el muestreo.
# Primero definimos la función plot_gauss, que nos permitirá graficar rápidamente.
plot_gauss <- function(x, name) {
  avg <- mean(x)
  std <- sd(x)
  x_axis <- seq(min(x), max(x), length = length(x))
  y_axis <- dnorm(x_axis, mean = avg, sd = std)
  hist(x, main = sprintf("Gráfico de la distribución %s", name))
  lines(x_axis, y_axis, type = "l", col = "red")
}

# Representamos las distribuciones Normal, Gamma y de Poisson, de parámetros N(1.5, 0.75), G(2, 1.5) y P(5).
avg <- 1.5
std <- 0.75
shape <- 2
scale <- 1.5
lambda <- 5

# Inicializamos los vectores que acumularán las medias para cada muestreo.
iter <- 1000
normal_vector <- vector(length = iter)
gamma_vector <- vector(length = iter)
poisson_vector <- vector(length = iter)

# Graficamos las 3 distribuciones para distintos tamaños de la muestra
sample_size <- 10

# Realizar el muestreo de las distribuciones.
for(i in 1:iter) {
  new_vector <- rnorm(n = sample_size, mean = avg, sd = std)
  normal_vector[i] <- mean(new_vector)
  
  new_vector <- rgamma(n = sample_size, shape = shape, scale = scale)
  gamma_vector[i] <- mean(new_vector)

  new_vector <- rpois(n = sample_size, lambda = lambda)
  poisson_vector[i] <- mean(new_vector)
}

# Código para gráficos
par(mfrow = c(1, 3))
plot_gauss(normal_vector, sprintf("Normal con tamaño %s", sample_size))
plot_gauss(gamma_vector, sprintf("Gamma con tamaño %s", sample_size))
plot_gauss(poisson_vector, sprintf("de Poisson con tamaño %s", sample_size))

# Graficamos las 3 distribuciones para distintos tamaños de la muestra
sample_size <- 100

# Realizar el muestreo de las distribuciones.
for(i in 1:iter) {
  new_vector <- rnorm(n = sample_size, mean = avg, sd = std)
  normal_vector[i] <- mean(new_vector)
  
  new_vector <- rgamma(n = sample_size, shape = shape, scale = scale)
  gamma_vector[i] <- mean(new_vector)

  new_vector <- rpois(n = sample_size, lambda = lambda)
  poisson_vector[i] <- mean(new_vector)
}

# Código para gráficos
par(mfrow = c(1, 3))
plot_gauss(normal_vector, sprintf("Normal con tamaño %s", sample_size))
plot_gauss(gamma_vector, sprintf("Gamma con tamaño %s", sample_size))
plot_gauss(poisson_vector, sprintf("de Poisson con tamaño %s", sample_size))

# Graficamos las 3 distribuciones para distintos tamaños de la muestra
sample_size <- 1000

# Realizar el muestreo de las distribuciones.
for(i in 1:iter) {
  new_vector <- rnorm(n = sample_size, mean = avg, sd = std)
  normal_vector[i] <- mean(new_vector)
  
  new_vector <- rgamma(n = sample_size, shape = shape, scale = scale)
  gamma_vector[i] <- mean(new_vector)

  new_vector <- rpois(n = sample_size, lambda = lambda)
  poisson_vector[i] <- mean(new_vector)
}

# Código para gráficos
par(mfrow = c(1, 3))
plot_gauss(normal_vector, sprintf("Normal con tamaño %s", sample_size))
plot_gauss(gamma_vector, sprintf("Gamma con tamaño %s", sample_size))
plot_gauss(poisson_vector, sprintf("de Poisson con tamaño %s", sample_size))

# Graficamos las 3 distribuciones para distintos tamaños de la muestra
sample_size <- 5000

# Realizar el muestreo de las distribuciones.
for(i in 1:iter) {
  new_vector <- rnorm(n = sample_size, mean = avg, sd = std)
  normal_vector[i] <- mean(new_vector)
  
  new_vector <- rgamma(n = sample_size, shape = shape, scale = scale)
  gamma_vector[i] <- mean(new_vector)

  new_vector <- rpois(n = sample_size, lambda = lambda)
  poisson_vector[i] <- mean(new_vector)
}

# Código para gráficos
par(mfrow = c(1, 3))
plot_gauss(normal_vector, sprintf("Normal con tamaño %s", sample_size))
plot_gauss(gamma_vector, sprintf("Gamma con tamaño %s", sample_size))
plot_gauss(poisson_vector, sprintf("de Poisson con tamaño %s", sample_size))

Notemos que para la distribución de Poisson, los datos convergen a una distribución Normal al realizar muestreos de tamaño 100, mientras que las otras dos aún no muestran ese comportamiento. Es al tamaño de muestra de 1000 que las 3 distribuciones se comportan como uno esperaría teóricamente.

Pregunta 3: Ley de los Grandes Numeros.

Lanzamiento de monedas

Realice el experimento de lanzar una moneda cargada 1000 veces y observe el comportamiento que tiene la probabilidad de salir cara. Para realizar el experimento considere que la moneda tiene una probabilidad de \(5/8\) de salir cara. Grafique el experimento para las secuencias de intentos que van desde 1 a 1000, señalando el valor en que converge la probabilidad de salir cara.

Respuesta

# Inicializa un vector para almacenar las probabilidades en cada paso
probabilidades <- numeric(1000)

# Simular lanzamientos
for (lanzamientos in 1:1000) {
  # Genera una secuencia de lanzamientos de la moneda y cuenta cuántas caras salieron
  resultados <- sample(c("cara", "sello"), lanzamientos, replace = TRUE, prob = c(5/8, 1 - 5/8))
  caras <- sum(resultados == "cara")
  
  # Calcula la probabilidad de salir cara en este punto
  probabilidad_actual <- caras / lanzamientos
  probabilidades[lanzamientos] <- probabilidad_actual
}

# Gráfico de la convergencia

plot(1:1000, probabilidades, type = "l", col = "blue", xlab = "Número de intentos", ylab = "Probabilidad de salir cara")
abline(h = 5/8, col = "red", lty = 2)
legend("topright", legend = c("Probabilidad teórica", "Probabilidad simulada"), col = c("red", "blue"), lty = 2:1)

El problema de Monty Hall

Remontándonos en la televisión del año 1963, en USA existía un programa de concursos donde los participantes debían escoger entre 3 puertas para ganar un premio soñado. El problema del concurso era que solo detrás de 1 puerta estaba el premio mayor, mientras que detrás de las otras dos habían cabras como “premio”.

Una de las particularidades de este concurso, es que cuando el participante escogía una puerta, el animador del show abría una de las puertas que no fue escogida por el participante (Obviamente la puerta abierta por el animador no contenía el premio). Tras abrir la puerta, el animador consultaba al participante si su elección era definitiva, o si deseaba cambiar la puerta escogida por la otra puerta cerrada. Un vídeo que puede ayudar a comprender mejor este problema en el siguiente link.

Imagine que usted es participante del concurso y desea calcular la probabilidad de ganar el gran premio si cambia de puerta en el momento que el animador se lo ofrece. Utilizando listas/arrays/vectores simule las puertas del concurso, dejando aleatoriamente el premio en alguna posición del array. Hecho esto, genere un numero de forma aleatoria para escoger una de las puerta (posiciones de la estructura), para luego ver si cambiando de posición tendrá mayores posibilidades de ganar el premio. Genere N veces el experimento y grafique cada una de las iteraciones, tal como se hizo en el ejercicio de las monedas.

Respuesta:

# Creamos una función que simule el juego
montyhall <- function(cambiar = TRUE){
  puertas <- sample(1:3,3)             #Puertas donde la posición que tiene el 3 es el premio
  eleccion <- sample(1:3,1)            #Elección del participante.

  puertas_no_elegidas <- setdiff(puertas, eleccion)

  puertas_no_elegidas_corr <-puertas_no_elegidas[puertas_no_elegidas != 3]
  if(length(puertas_no_elegidas_corr) == 1){
    puerta_abierta <- puertas_no_elegidas[1]
  }
  else{
    puerta_abierta <- sample(puertas_no_elegidas[puertas_no_elegidas != 3], size = 1)
  }
  


  eleccion_final <- eleccion == 3
  if(cambiar){

    puertas_restantes <- setdiff(puertas, c(eleccion, puerta_abierta))

    eleccion_final <- puertas_restantes
    eleccion_final <- puertas[eleccion_final] == 3
  }
  
  return(eleccion_final) # Retornamos la elección, esta puede que tenga el premio o no
}

# Función que simula N juegos
n_juegos <- function(n = 10 ,cambiar_puerta = TRUE){
  resultados <- replicate(n, montyhall(cambiar_puerta))
  
  probabilidad_cambio <- mean(resultados)
  
  plot(1:n, cumsum(resultados)/(1:n), type = "l", col = "blue", 
     xlab = "Número de simulaciones", ylab = "Probabilidad de ganar ")
abline(h = probabilidad_cambio, col = "red", lty = 2)
legend("topright", legend = c("Probabilidad", "Probabilidad promedio"), col = c("red", "blue"), lty = 2:1)
  

}
n_juegos(1000)

n_juegos(1000,FALSE)

Pregunta 4: ¿Independencia?

Ustedes disponen de los dados D1 y D2, los cuales no están cargados y son utilizados para comprobar que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\) cuando el evento A es independiente del B. Para estudiar la independencia considere que los eventos A y B se definen de la siguiente manera; sea A el evento dado por los valores obtenidos en el lanzamiento del dado D1, este está compuesto por \(A=\{D1=1,D1=2,D1=6\}\). Por otro lado, el evento B viene dado por los valores obtenidos con el dado D2, el que está conformado por \(B=\{D2=1,D2=2,D2=3,D2=4\}\). Con esto, tendremos un \(\mathbb{P}(A)=1/2\), \(\mathbb{P}(𝐵)=2/3\) y \(\mathbb{P}(AB)=1/3\). Compruebe de forma gráfica que al realizar 1000 lanzamientos (u otro valor grande que usted desea probar) se visualiza que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\).

Hecho lo anterior, compruebe el comportamiento de un segundo grupo de eventos, dados por el lanzamiento de solo el dado D1. Donde, los eventos para D1 quedan definidos como: \(A =\{D1=1,D1=2,D1=6\}\) y \(B=\{D1=1,D1=2,D1=3\}\). ¿Se observa independencia en este experimento?.

Se le aconseja que para simular los lanzamientos de dados utilice la función sample() para generar valores aleatorios entre 1 y 6. Compruebe los números generados por la función con los casos favorables de cada uno de los eventos a ser estudiados.

Respuesta:

n_lan <- 6000000

# Primer grupo de eventos
dice1 <- sample(1:6, n_lan, replace = TRUE)  # Dado 1
dice2 <- sample(1:6, n_lan, replace = TRUE)  # Dado 2
dice_count1 <- table(dice1)
dice_count2 <- table(dice2)

favorable_a <- c(1, 2, 6)  # {1, 2, 6}
favorable_b <- c(1, 2, 3, 4)  # {1, 2, 3, 4}
favorable_ayb <- c(1, 2)  # {1, 2}

count_a <- sum(dice_count1[favorable_a])  # |A|
count_b <- sum(dice_count2[favorable_b])  # |B|

# Ahora contamos los casos donde A y B se cumplen al mismo tiempo
count_ayb <- 0  # |A∩B|
for (i in 1:n_lan) {
  if (dice1[i] %in% favorable_a && dice2[i] %in% favorable_b) {
    count_ayb <- count_ayb + 1
  }
}

# Calculamos las probabilidades
prob_a <- count_a / n_lan  # P(A)
prob_b <- count_b / n_lan  # P(B)
prob_ayb <- count_ayb / n_lan  # P(A∩B)

print(sprintf("P(A) = %s", prob_a))

## [1] "P(A) = 0.4999925"

print(sprintf("P(B) = %s", prob_b))

## [1] "P(B) = 0.666550833333333"

print(sprintf("P(A)*P(B) = %s", prob_a*prob_b))

## [1] "P(A)*P(B) = 0.333270417535417"

print(sprintf("P(A∩B) = %s", prob_ayb))

## [1] "P(A∩B) = 0.333211833333333"

Ahora visualizamos las probabilidades

par(mfrow = c(1, 2))
y_ax <- c(0, 0, 1, 1)
x_ax <- c(0, 1, 1, 0)

# Graficamos P(A)
plot(1, ann = FALSE, type = "n", axes = FALSE, xlim = c(0, 1), ylim = c(0, 1), main = "P(A)")
x_p_a <- c(0, prob_a, prob_a, 0)
x_p_a_c <- c(prob_a, 1, 1, prob_a)

polygon(x_p_a, y_ax, col = "blue")
text(mean(x_p_a), mean(y_ax), expression(P(A)), col = "white", cex = 1.5)

polygon(x_p_a_c, y_ax, col = "red")
text(mean(x_p_a_c), mean(y_ax), expression(P(A^c)), col = "white", cex = 1.5)


# Graficamos P(B)
plot(2, ann = FALSE, type = "n", axes = FALSE, xlim = c(0, 1), ylim = c(0, 1), main = "P(B)")
y_p_b <- c(prob_b, prob_b, 1, 1)
y_p_b_c <- c(0, 0, prob_b, prob_b)

polygon(x_ax, y_p_b, col = "blue")
text(mean(x_ax), mean(y_p_b), expression(P(B)), col = "white", cex = 1.5)

polygon(x_ax, y_p_b_c, col = "red")
text(mean(x_ax), mean(y_p_b_c), expression(P(B^c)), col = "white", cex = 1.5)

Ahora grafiquemos \(\mathbb{P}(AB)\)

plot(1, ann = FALSE, type = "n", axes = FALSE, xlim = c(0, 1), ylim = c(0, 1), main = "P(A∩B)")

polygon(c(0, prob_a, prob_a, 0), c(prob_b, prob_b, 1, 1), col = "blue")
text(prob_a / 2, 1 - prob_b / 4, expression(P(A * "∩" * B)), col = "white", cex = 1.5)

polygon(c(0, prob_a, prob_a, 0), c(0, 0, prob_b, prob_b), col = "purple")
text(prob_a / 2, prob_b / 2, expression(P(A * "∩" * B^c)), col = "white", cex = 1.5)

polygon(c(prob_a, 1, 1, prob_a), c(prob_b, prob_b, 1, 1), col = "purple")
text(1 - prob_a / 2, 1 - prob_b / 4, expression(P(A^c * "∩" * B)), col = "white", cex = 1.5)

polygon(c(prob_a, 1, 1, prob_a), c(0, 0, prob_b, prob_b), col = "red")
text(1 - prob_a / 2, prob_b / 2, expression(P(A^c * "∩" * B^c)), col = "white", cex = 1.5)

Notemos que el área del cuadrado de la esquina superior izquierda, que representa el área de \(\mathbb{P}(A \cap B)\), tiene la misma área que la superposición de las figuras anteriores,

Pregunta 5: La Ruina del Jugador

Un amigo ludópata suyo le comenta que el truco de jugar en el casino esta en no parar de apostar y apostando lo mínimo posible. Ya que así, tienes mas probabilidades de ganar el gran pozo que acumula el juego. Usted sabiendo la condición de su amigo, decide no creer en su conjetura y decide probar esto a través de un experimento.

Para realizar el experimento usted decide asumir las siguientes declaraciones, bajo sus observaciones:

La probabilidad de ganar en un juego del casino es \(8/19\)
Sabe que su amigo solo posee siempre fondos en el rango de 0 a 200 dolares.
Las apuestas como mínimo deben ser igual a 5 dolares.
El monto de las apuestas no cambia y son siempre igual a la primera. Por ejemplo, si su amigo apuesta 50 dolares, todos los próximos juegos apuesta 50 hasta que se acaba su dinero.
Usted pierde cuando se le acaban los fondos.

En el primer experimento deberá obtener la evolución de los fondos hasta que el jugador se queda sin fondos para jugar. Puede ser útil seguir la lógica de una moneda cargada para realizar esto. Pruebe esto con una apuesta igual a 5, 25 y 50 graficando los resultados. Comente los resultados obtenidos.

Para la segunda parte de este experimento, con las funciones ya creadas, proyecte 5000 apuestas y obtenga la probabilidad a la que converge el experimento empezando con una apuesta de: 5, 25 y 50. Para este experimento considerar como éxito (1) cuando su función ruina supera los 200 y considere lo contrario cuando su función perdida es menor o igual a 0.

Respuesta

# Función para obtener el desarrollo de las apuestas
ruina <- function(money = 100, bet = 5) {
  vec_money <- c(money)
  while (0 < money && money < 200) {
    game <- as.integer(runif(1, 1, 19))
    sign <- if (game <= 8) 1 else -1
    money <- money + sign * bet
    vec_money <- append(vec_money, money)
  }
  return(vec_money)
}

plot(ruina(), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 5)")

plot(ruina(bet = 25), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 25)")

plot(ruina(bet = 50), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 50)")

En general se ve que la tendencia es perder, solo que al apostar el mínimo se pierde más lento ya que al ser una variación más pequeña, se desestabiliza (es decir, se va a los extremos) de a poco, no así al apostar 50, que es posible perder en dos juegos.

n_bet <- 5000
bet5_vector <- vector(length = n_bet)
bet20_vector <- vector(length = n_bet)
bet50_vector <- vector(length = n_bet)
for (i in 1:n_bet) {
  final5_money <- tail(ruina(bet = 5), 1)
  final20_money <- tail(ruina(bet = 20), 1)
  final50_money <- tail(ruina(bet = 50), 1)
  win5 <- if (final5_money <= 0) 0 else 1
  win20 <- if (final20_money <= 0) 0 else 1
  win50 <- if (final50_money <= 0) 0 else 1
  bet5_vector[i] <- win5
  bet20_vector[i] <- win20
  bet50_vector[i] <- win50
}
t5 <- table(bet5_vector)
t20 <- table(bet20_vector)
t50 <- table(bet50_vector)
print(sprintf("Apuesta = 5 -> E = %d, F = %d, Prob = %f", t5["1"], t5["0"], t5["1"]/(t5["1"] + t5["0"])))

## [1] "Apuesta = 5 -> E = 54, F = 4946, Prob = 0.010800"

hist(bet5_vector, main = sprintf("Éxitos y fracasos para apuesta = 5"))

print(sprintf("Apuesta = 20 -> E = %d, F = %d, Prob = %f", t20["1"], t20["0"], t20["1"]/(t20["1"] + t20["0"])))

## [1] "Apuesta = 20 -> E = 1217, F = 3783, Prob = 0.243400"

hist(bet20_vector, main = sprintf("Éxitos y fracasos para apuesta = 20"))

print(sprintf("Apuesta = 50 -> E = %d, F = %d, Prob = %f", t50["1"], t50["0"], t50["1"]/(t50["1"] + t50["0"])))

## [1] "Apuesta = 50 -> E = 1959, F = 3041, Prob = 0.391800"

hist(bet50_vector, main = sprintf("Éxitos y fracasos para apuesta = 50"))

Puede verse que la probabilidad de éxito va aumentando a medida que el monto de apuesta aumenta, ya que es más fácil ganar al ser menor la cantidad de juegos necesarios para pasar la barrera de los 200.

A work by CC6104